Motion artefact correction and physiological noise removal

This work is based on the excellent tutorials created by the MNE and MNE-NIRS developers:

• MNE: Tutorials
• MNE-NIRS: Examples

Table of Contents

• Pre-processing
  • Load python modules
  • Import data
    • Import snirf data
    • Add info to raw intensity data
  • Signal quality metrics and channel rejection
    • Visual inspection
    • Scalp-coupling index
    • Peak spectral power
  • Motion artefact correction
  • Short channel regression
    • Method of Saager and Berger
  • Global component derived from the mean or median
    • Method of Saager and Berger
  • Signal processing methods
    • Lowpass filtering

Pre-processing

---

Load python modules

# Import the needed modules 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import mne
import mne_nirs
import pywt
import scipy
import dabest
import pathlib
from itertools import compress
from mne.preprocessing.nirs import scalp_coupling_index, optical_density, temporal_derivative_distribution_repair, beer_lambert_law, _validate_nirs_info, source_detector_distances
from mne_nirs.preprocessing import peak_power, scalp_coupling_index_windowed, quantify_mayer_fooof
from mne_nirs.visualisation import plot_timechannel_quality_metric
from mne_nirs.channels import get_short_channels, get_long_channels, picks_pair_to_idx
from mne_nirs.signal_enhancement import short_channel_regression
from mne_nirs.experimental_design import make_first_level_design_matrix
from mne.preprocessing import ICA, corrmap
from mne_bids import BIDSPath, read_raw_bids
from scipy.stats import norm, ttest_1samp, kurtosis, pearsonr
from scipy.signal import firwin, freqz, filtfilt
from scipy.interpolate import CubicSpline
from fooof import FOOOF
from csaps import csaps
from sklearn.decomposition import PCA, FastICA
from sklearn.preprocessing import StandardScaler

C:\ProgramData\anaconda3\Lib\site-packages\paramiko\transport.py:219: CryptographyDeprecationWarning: Blowfish has been deprecated and will be removed in a future release
  "class": algorithms.Blowfish,
C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\1896447934.py:24: DeprecationWarning: 
The `fooof` package is being deprecated and replaced by the `specparam` (spectral parameterization) package.
This version of `fooof` (1.1) is fully functional, but will not be further updated.
New projects are recommended to update to using `specparam` (see Changelog for details).
  from fooof import FOOOF

# Write own fucntions
def reorder(raw_initial):
    """
    Input: raw object
    Function: reorders the channels to match the MNE standard
    """
    raw = raw_initial.copy()
    channel_names = raw.ch_names
    channel_names_ordered = []
    for i in range(0, int(len(channel_names)/2)):
        channel_names_ordered.append(channel_names[i])
        channel_names_ordered.append(channel_names[(i + int(len(channel_names)/2))])
    raw_ordered = raw.reorder_channels(channel_names_ordered)
    return raw_ordered

def add_info(raw_intensity, first_name, last_name, age, sex, EHI, annotations, duration_annotation, duration_rest, N_control):
    """
    Input: raw object, first name, last name, age, sex, EHI, annotations, duration of annotation, duration of baseline period, # of control trials
    Function: Returns raw object with added subject info and random allocated control trials within the baseline period
    """
    raw = raw_intensity.copy()
    
    # Add subject info
    raw.info["subject_info"] = {'first_name' : first_name, 'last_name' : last_name, 'age' : age, 'sex' : sex, 'EHI' : EHI}
    
    # attribute meaningful names to trigger codes. Include information about the duration of each stimulus (10 seconds). 
    raw.annotations.set_durations(duration_annotation)
    raw.annotations.rename(annotations)
    
    events, event_dict = mne.events_from_annotations(raw)
    
    # Define begin and end of rest period
    begin_rest = np.ceil(events[np.where(events[:,2] == event_dict['Baseline'])[0][0],0]/raw.info['sfreq'])
    end_rest = begin_rest + duration_rest - duration_annotation 
    
    # Remove start indicator of rest period
    raw.annotations.delete(np.nonzero(raw.annotations.description == "Baseline"))
    
    # Define indices of random selected control trials
    indices = np.random.choice(np.arange(begin_rest, end_rest), N_control, replace = False).astype(int)
    raw.annotations.append(indices, [duration_annotation]*N_control, ["Baseline"]*N_control)
    
    return raw

def add_info2(raw_intensity, first_name, last_name, age, sex, EHI):
    """
    Input: raw object, first name, last name, age, sex, EHI
    Function: Returns raw object with added subject info
    """
    raw = raw_intensity.copy()
    
    # Add subject info
    raw.info["subject_info"] = {'first_name' : first_name, 'last_name' : last_name, 'age' : age, 'sex' : sex, 'EHI' : EHI}
    
    return raw
    
def crop(raw_input, t_before = 10, t_after = 60):
    """
    Input: raw object, time interval length before first event and time interval length after last event that should be included in the recording
    Function: crops raw object to a recording starting 't_before' s (default 10 s) before the first annotated event and ending 't_after' s (default 60 s) after the last event 
    """
    raw = raw_input.copy()
    end_rec = raw.get_data().shape[1]/raw.info['sfreq']
    t_start = events[0,0]/raw.info['sfreq'] - t_before
    t_end = events[-1,0]/raw.info['sfreq'] + t_after
    if t_start < 10:
        t_start = 0
    if t_end > end_rec:
        t_end = end_rec
    raw.crop(t_start, t_end)
    return raw

def get_indices(array, target_array):
    """
    Input: complete array of channel names, array of channels names to be removed
    Function: gives back indices of desired channels
    """
    lst = list(array)
    indices = list(np.arange(0, len(array)))
    for target in target_array:
        if target in lst:
            indices.remove(lst.index(target))
    return indices

def ch_names_L_R(raw):
    """
    Input: raw object
    Function: Separate list of channel names in 2 lists for the left and right hemisphere respectively
    """
    left, right = [], []
    for ch in raw.ch_names:
        if (int(ch[1]) % 2) == 0:
            right.append(ch)
        else:
            left.append(ch)
    return left, right

def scale_up_spectra(spectra, freqs):
    """
    Input: spectra, freqs
    Function: FOOOF requires the frequency values to be higher than the fNIRS data permits, 
    so we scale the values up by 10 here, and then will scale
    the frequency values down by 10 later.
    """
    freqs = freqs * 10
    return spectra, freqs

def get_IMU_data(filename, plot = False):
    """
    Extract IMU data from text file.
    Input: filename with .txt, plot (bool)
    Output: Three dataframes containing the complete IMU data, the accelerometer data and the gyroscope data respectively
    """
    
    # Read .txt file and convert to dataframe
    df = pd.read_csv(filename, sep=';', header=None, names=["time", "battery", "channels", "gyroX", "gyroY", "gyroZ", "accX", "accY", "accZ", "marker", "_"])

    # Select IMU data
    data = df.drop(0)
    markers = data["marker"].astype(float)
    IMU_data = data.drop(columns =['time', 'battery', 'channels', 'marker', '_']).astype(float)

    gyro_data = IMU_data.drop(columns = ['accX', 'accY', 'accZ'])
    acc_data = IMU_data.drop(columns = ["gyroX", "gyroY", "gyroZ"])
    
    if plot: 
        %matplotlib inline
        gyro_data.plot(title='Gyroscope data', grid = True, xlabel = 'samples', ylabel = 'dps')
        acc_data.plot(title='Accelerometer data', grid = True, xlabel = 'samples', ylabel = 'm/s^2')

    return IMU_data, acc_data, gyro_data

def compare_original_filtered(raw_filtered, raw_unfiltered, export = False, filename = None):
    data_new = raw_filtered.get_data()
    data = raw_unfiltered.get_data()
    picks = np.sort(_validate_nirs_info(raw_filtered.info))
    datatype = raw_filtered.ch_names[0][-3:]
    for pick in picks:
        if plot:
            if datatype == 'hbo':
                type = ' hbo' if pick%2 == 0 else ' hbr'
                ylabel = 'HbO/HbR (M)'
            else:
                type = ' wavelength 760' if pick%2 == 0 else ' wavelength 850'
                ylabel = 'OD (V)'
            %matplotlib inline
            time = np.arange(data.shape[1])/raw_filtered.info['sfreq']
            fig = plt.figure(figsize = (12,5))
            plt.plot(time, data[pick], label = 'original')
            plt.plot(time, data_new[pick], label = 'filtered')
            plt.xlabel('time (s)')
            plt.ylabel(ylabel)
            plt.title('Channel: ' + str(pick//2+1) + type)
            plt.legend()
            if export:
                if filename == None:
                    raise ValueError('Filename must be given to export figures')
                fig.savefig(filename + '_channel_' + str(pick//2+1) + type + '.png')
            plt.show()

# Set plot characteristics
%matplotlib inline
sns.set_theme() # nicer plots
# Load external plots?
plot = False
# Export figures?
export = False

Import data

Import snirf data

raws = []
data_dir = pathlib.Path("C:/Users/fabia/fNIRS data analysis/N-back analysis/N_back data_BIDS")
bids_root = data_dir.with_name(data_dir.name)
for sub in range(1, 5):  # Loop from first to fourth subject
    # Create path to file based on experiment info
    bids_path = BIDSPath(subject="%02d" % sub, task="Nback", datatype="nirs",
                         root=bids_root, suffix="nirs",
                         extension=".snirf")
    raws.append(read_raw_bids(bids_path=bids_path, verbose=False))
raws

[<RawSNIRF | sub-01_task-Nback_nirs.snirf, 32 x 20463 (2046.2 s), ~36 kB, data not loaded>,
 <RawSNIRF | sub-02_task-Nback_nirs.snirf, 32 x 40258 (4549.0 s), ~36 kB, data not loaded>,
 <RawSNIRF | sub-03_task-Nback_nirs.snirf, 32 x 20339 (2277.9 s), ~36 kB, data not loaded>,
 <RawSNIRF | sub-04_task-Nback_nirs.snirf, 32 x 21552 (2650.8 s), ~36 kB, data not loaded>]

raw_intensities = []
# Reorder channels (Due to non-uniformity between mne_nirs and .snirf data)
for sub in range(4):
    raw_intensities.append(reorder(raws[sub]))

Add info to raw intensity data

raw_intensities[0] = add_info2(raw_intensities[0], 'Subject', '1', 23, 'male', 'right')
raw_intensities[1] = add_info2(raw_intensities[1], 'Subject', '2', 22, 'female', 'right')
raw_intensities[2] = add_info2(raw_intensities[2], 'Subject', '3', 25, 'male', 'right')
raw_intensities[3] = add_info2(raw_intensities[3], 'Subject', '4', 25, 'female', 'right')

# Look at characteristics
sfreqs = []
lowpass = []
for sub in range(4):
    sfreqs.append(raw_intensities[sub].info['sfreq'])
    lowpass.append(raw_intensities[sub].info['lowpass'])
print('Sample frequencies (Hz): ' + str(sfreqs))
print('Max bandwidth (Hz): ' + str(lowpass))

Sample frequencies (Hz): [10.0, 8.849557522123893, 8.928571428571429, 8.130081300813009]
Max bandwidth (Hz): [5.0, 4.424778761061947, 4.464285714285714, 4.065040650406504]

# Look at segments
events = []
event_dicts = []
for sub in range(4):
    event, event_dict = mne.events_from_annotations(raw_intensities[sub])
    events.append(event)
    event_dicts.append(event_dict)
%matplotlib inline
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(21, 11))

mne.viz.plot_events(events[0], event_id=event_dicts[0], sfreq=raw_intensities[0].info['sfreq'], axes=axes[0,0], show = False)
mne.viz.plot_events(events[1], event_id=event_dicts[1], sfreq=raw_intensities[1].info['sfreq'], axes=axes[0,1], show = False)
mne.viz.plot_events(events[2], event_id=event_dicts[2], sfreq=raw_intensities[2].info['sfreq'], axes=axes[1,0], show = False)
mne.viz.plot_events(events[3], event_id=event_dicts[3], sfreq=raw_intensities[3].info['sfreq'], axes=axes[1,1], show = False)
axes[0,0].set_title('Participant 1', fontweight="bold")
axes[0,1].set_title('Participant 2', fontweight="bold")
axes[1,0].set_title('Participant 3', fontweight="bold")
axes[1,1].set_title('Participant 4', fontweight="bold")

fig.show()

Used Annotations descriptions: ['0Back', '2Back_block1', '2Back_block2', '2Back_block3', 'Baseline', 'Practice/0Back', 'Practice/2Back']
Used Annotations descriptions: ['0Back', '2Back_block1', '2Back_block2', '2Back_block3', 'Baseline', 'Practice/0Back', 'Practice/2Back']
Used Annotations descriptions: ['0Back', '2Back_block1', '2Back_block2', '2Back_block3', 'Baseline', 'Practice/0Back', 'Practice/2Back']
Used Annotations descriptions: ['0Back', '2Back_block1', '2Back_block2', '2Back_block3', 'Baseline', 'Practice/0Back', 'Practice/2Back']

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\4025079886.py:20: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()

---

Remark:

• 16 channels * 2 wavelengths = 32 channels
• Different sample frequencies are used for each participant = [10.0, 8.849557522123893, 8.928571428571429, 8.130081300813009] Hz
• LPF = [5.0, 4.424778761061947, 4.464285714285714, 4.065040650406504] Hz - Related to max bandwidth (Nyquist theorem): $F_{max} = \frac{F_s}{2}$
• 10 practice 0-back and 2-back letters were presented
• 40 letters were presented for the 0-back and the baseline recording
• 400 letters were presented for each of the three 2-back blocks

---

Validate that the location of sources-detector pairs and channels are in the expected locations

# Plot channels in 2D using build-in function of mne_nirs
sns.reset_defaults() # Turn of sns to be able to also plot contour of head
# dots = Center of channels
%matplotlib inline
raw_intensities[0].plot_sensors();

sns.set_theme()

---

Remark:

• 14 long channels and 2 short channels
• 12 sources, 4 detectors and 2 short channel sources
• montage located above prefrontal cortex

---

Signal quality metrics and channel rejection

Visual inspection

Visual inspection of the raw light intensity data

Subject 1

# Plot raw data with build-in plot function of mne_nirs
if plot:
    %matplotlib qt
    raw_intensities[0].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 2

# Plot raw data with build-in plot function of mne_nirs
if plot:
    %matplotlib qt
    raw_intensities[1].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 3

# Plot raw data with build-in plot function of mne_nirs
if plot:
    %matplotlib qt
    raw_intensities[2].plot(duration=300, show_scrollbars = True, clipping = None);

---

Remark:

• Channel S14-D4 and S12-D2 seem to contain no data at all. Important as S14-D4 is a short channel.

---

Subject 4

# Plot raw data with build-in plot function of mne_nirs
if plot:
    %matplotlib qt
    raw_intensities[3].plot(duration=300, show_scrollbars = True, clipping = None);

Visual inspection of the the optical density

The OD or the attenuation of incident light, can be calculated as the logarithmic ratio of the input light intensity ($I_{in}$) and the detected light intensity ($I_{out}$): $$OD_\lambda = \log \left(\frac{I_{in}}{I_{out}}\right)$$

# Convert raw intensity data to OD with build-in function of mne_nirs
raw_ods = []
for sub in range(4):
    raw_ods.append(optical_density(raw_intensities[sub]))

Reading 0 ... 20462  =      0.000 ...  2046.200 secs...
Reading 0 ... 40257  =      0.000 ...  4549.041 secs...
Reading 0 ... 20338  =      0.000 ...  2277.856 secs...
Reading 0 ... 21551  =      0.000 ...  2650.773 secs...

Remark: Negative values in raw intensity due to artefacts

Subject 1

# Plot OD
if plot:
    %matplotlib qt
    raw_ods[0].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 2

# Plot OD
if plot:
    %matplotlib qt
    raw_ods[1].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 3

# Plot OD
if plot:
    %matplotlib qt
    raw_ods[2].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 4

# Plot OD
if plot:
    %matplotlib qt
    raw_ods[3].plot(duration=300, show_scrollbars = True, clipping = None);

---

Remark:

• For more information about the visual inspection procedure can be found in the Signal quality metrics and channel rejection - N-back Notebook.

---

Separate long and short channels

The short channels were not assessed using signal quality metrics because they inherently contain physiological noise. These channels cannot be removed from the dataset as they are essential for performing short channel regression.

raw_ods_full = []
raw_ods_short = []
raw_ods_new = []
min_dist = 0.01 
max_dist = 0.045
for sub in range(4):
    raw_ods_full.append(raw_ods[sub].copy())
    raw_ods_short.append(get_short_channels(raw_ods[sub], max_dist=min_dist))
    raw_ods_new.append(get_long_channels(raw_ods[sub], min_dist=min_dist, max_dist=max_dist))

raw_ods = raw_ods_new

# Short channels are: (channel number 11 and 12)
raw_ods_full[0].ch_names[20:24]

['S13_D3 760', 'S13_D3 850', 'S14_D4 760', 'S14_D4 850']

Scalp-coupling index

The SCI is a quality measurement of the connection between the optodes and the scalp of the participant using the strong cardiac oscillation in raw fNIRS signals. With a sufficiently high sampling rate (e.g.: 10 Hz), the heartbeat emerges as a dependable marker for evaluating optode-scalp coupling. It was first introduced by Pollonini et al. (2013) and it assesses the synchronous cardiac pulsation in the two photodetected signals (760 and 850 nm) in the cardiac band (0.5 - 2.5 Hz). After the bandpass-filtering the resulting signals are normalized to balance any difference between their amplitude. Finally, the zero-lag cross-correlation between photodetected signals of the same channel is used as a quantitative measure of the SNR of that specific channel. A threshold of SCI $>$ 0.75 is recommended by the authors.

Given that the frequency band of 0.5 - 2.5 Hz is quite broad and the participants are relatively at rest while seated, a FIR filter with a narrower band of 0.7 - 1.5 Hz is applied (40 bpm - 85 bpm). The transition band width at both the low and high cut-off frequencies is set to 0.3 Hz.

SCI evaluated over complete signal

# Calculate SCI with build-in function of mne_nirs
scis = []
for sub in range(4):
    scis.append(scalp_coupling_index(raw_ods[sub], l_freq=0.7, h_freq=1.5, l_trans_bandwidth=0.3, h_trans_bandwidth=0.3, verbose=False))
    # SCI bandpasses the signal with a FIR filter: l_freq=0.7, h_freq=1.5, l_trans_bandwidth=0.3, h_trans_bandwidth=0.3

# Copy raw object to apply SCI 
raw_od_SCIs = []
for sub in range(4):
    raw_od_SCIs.append(raw_ods[sub].copy())
# Mark channels with SCI < 0.75 as BAD
for sub in range(4):
    raw_od_SCIs[sub].info['bads'] = list(compress(raw_ods[sub].ch_names, scis[sub] < 0.75)) # 0.75 RECOMMENDED BY Pollonini et al. (2013)
    print('The ' + str(len(raw_od_SCIs[sub].info['bads'])) + ' bad channels for subject ' + str(sub+1) + ' are: ' + str(raw_od_SCIs[sub].info['bads']))

The 10 bad channels for subject 1 are: ['S4_D3 760', 'S4_D3 850', 'S6_D3 760', 'S6_D3 850', 'S7_D4 760', 'S7_D4 850', 'S8_D3 760', 'S8_D3 850', 'S9_D4 760', 'S9_D4 850']
The 12 bad channels for subject 2 are: ['S4_D3 760', 'S4_D3 850', 'S5_D4 760', 'S5_D4 850', 'S6_D3 760', 'S6_D3 850', 'S7_D4 760', 'S7_D4 850', 'S8_D3 760', 'S8_D3 850', 'S9_D4 760', 'S9_D4 850']
The 4 bad channels for subject 3 are: ['S9_D4 760', 'S9_D4 850', 'S12_D2 760', 'S12_D2 850']
The 18 bad channels for subject 4 are: ['S2_D2 760', 'S2_D2 850', 'S5_D4 760', 'S5_D4 850', 'S6_D3 760', 'S6_D3 850', 'S7_D4 760', 'S7_D4 850', 'S8_D3 760', 'S8_D3 850', 'S9_D4 760', 'S9_D4 850', 'S11_D2 760', 'S11_D2 850', 'S12_D1 760', 'S12_D1 850', 'S12_D2 760', 'S12_D2 850']

SCI evaluated over moving window

Some extra channels are marked bad due to consistent bad scores over time (see Signal quality metrics and channel rejection - N-back Notebook).

# Add bad channels
raw_od_SCIs[3].info['bads'] += ['S2_D2 760', 'S2_D2 850']

---

Remark:

• For more information about the scalp-oupling index procedure can be found in the Signal quality metrics and channel rejection - N-back Notebook.

---

Peak spectral power

Pollonini et al. (2016) state that to enhance the reliability of a quantitative measure for scalp coupling, the spectral power of the cross-correlated signal can serve as an additional indicator of cardiac signal strength. By setting a spectral power threshold, it becomes possible to objectively differentiate between a clear optical channel and a noisy one. Theoretically, photodetected cardiac signals can be conceptualized as two sinusoidal waves oscillating at the physiological frequency $f_{cardiac}$ (e.g.: 1 Hz or 60 bpm) with arbitrary amplitude and offset. The resultant normalized cross-correlation signal ($\overline{x_{\lambda_1}} \otimes \overline{x_{\lambda_2}}$) exhibits unity amplitude, a frequency of $f_{cardiac}$ and a peak power of 0.5. Therefore, a threshold can be established at a certain fraction of the ideal peak power. However, in practice, cardiac photoplethysmographic signals are not strictly sinusoidal and may contain quantization noise, which reduces the peak power value. Hence, the authors recommend empirically setting the threshold value at 0.1.

# View instances where a subset of channels may be contaminated by artifacts for a short duration of the recording
raw_od_SCI_PSPs = []
for sub in range(4):
    raw_od_SCI_PSP, scores, times = peak_power(raw_od_SCIs[sub], time_window=10)
    raw_od_SCI_PSPs.append(raw_od_SCI_PSP)

---

Remark:

• For more information about the peak spectral power procedure can be found in the Signal quality metrics and channel rejection - FTT 2 Notebook.

---

Check if signal quality assesment was performed well

Subject 1

# Plot OD
if plot:
    %matplotlib qt
    raw_od_SCI_PSPs[0].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 2

# Plot OD
if plot:
    %matplotlib qt
    raw_od_SCI_PSPs[1].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 3

# Plot OD
if plot:
    %matplotlib qt
    raw_od_SCI_PSPs[2].plot(duration=300, show_scrollbars = True, clipping = None);

Subject 4

# Plot OD
if plot:
    %matplotlib qt
    raw_od_SCI_PSPs[3].plot(duration=300, show_scrollbars = True, clipping = None);

---

Remark:

• Some time segments are marked manually as bad
• Selection of channels marked as bad through SCI
• Peak spectral power marked time segments in selected channels bad

---

Motion artefact correction

Temporal derivative distribution repair

The TDDR method was introduced by Fishburn et al. in 2019. It offers both online and offline filtering approaches based on the temporal derivative of fNIRS signals, requiring no user-defined parameters. Moreover, it is applicable to concentration changes, optical intensities, and optical densities (Huang et al., 2022).

This method operates under specific assumptions: (1) non-motion-related fluctuations follow a normal distribution, (2) the majority of fluctuations are unrelated to motion artifacts, and (3) MA derivatives dominate in the derivatives of fNIRS signals during their presence.

The algorithm comprises five distinct steps. Fishburn et al. (2019) provides an in-depth description of each step.

Perform TDDR on OD data

# Add short channels again
for sub in range(4):
    raw_od_SCIs[sub].add_channels([raw_ods_short[sub]])

raw_ods_TDDR = []
for sub in range(4):
    raw_ods_TDDR.append(temporal_derivative_distribution_repair(raw_od_SCIs[sub]))

Visualize results

# Show subject one as example for html page
plot = True

Subject 1

if plot:
    compare_original_filtered(raw_ods_TDDR[0], raw_ods_full[0], export = export, filename = 'TDDR/Subject_1')

plot = False

Subject 2

if plot:
    compare_original_filtered(raw_ods_TDDR[1], raw_ods_full[1], export = export, filename = 'TDDR/Subject_2')

Subject 3

if plot:
    compare_original_filtered(raw_ods_TDDR[2], raw_ods_full[2], export = export, filename = 'TDDR/Subject_3')

Subject 4

if plot:
    compare_original_filtered(raw_ods_TDDR[3], raw_ods_full[3], export = export, filename = 'TDDR/Subject_4')

raw_ods_TDDR[1].ch_names[20]

'S10_D1 760'

---

Remark:

• TDDR succeeds in removing most spikes
• Baseline drifts are still present
• Amplitude of some channels is changed

---

Short channel regression

Method of Saager and Berger

The reflectance measurements at two detectors, L and S, allow for the calculation of absorbance changes $\Delta A_L (t)$ and $\Delta A_S (t)$. To eliminate the influence of the superficial component of the probed tissue, measurements can be taken using two SDD: long and short. This can be achieved using the equation: $$ \Delta A_C (\lambda, t) = \Delta A_L (\lambda, t) - \alpha \cdot \Delta A_S (\lambda, t)$$

Here, $\Delta A_C (\lambda, t)$ represents the corrected attenuation, $\Delta A_L (\lambda, t)$ corresponds to the attenuation at the long SDD (sampling both brain and superficial components), and $\Delta A_S (\lambda, t)$ pertains to the short SDD (sampling only superficial components) (Scholkmann et al., 2014).

According to Saager and Berger (2005), the scaling parameter $\alpha$ can be determined directly using a least-squares method:

$$\alpha = \frac{\Delta A_S (\lambda, t) \cdot \Delta A_L (\lambda, t)}{\Delta A_S (\lambda, t) \cdot \Delta A_S (\lambda, t)}$$

# Short channels predominantly contain systemic responses and long channels have both neural and systemic contribution
raw_ods_SCR = []
for sub in range(4):
    raw_ods_SCR.append(short_channel_regression(raw_ods_TDDR[sub].copy(), max_dist=0.01))

C:\Users\fabia\AppData\Roaming\Python\Python311\site-packages\mne_nirs\signal_enhancement\_short_channel_correction.py:64: RuntimeWarning: invalid value encountered in double_scalars
  alfa = np.dot(A_s, A_l) / np.dot(A_s, A_s)

---

Remark:

• Error arises for subject 3 due to the fact that short channel S14-D4 contains no data (array filled with zeros). As a result, the dot-product of A_s in the SCR algorithm delivers zero values, which leads to division by zero.

---

Visualize result

plot = True

Subject 1

if plot:
    compare_original_filtered(raw_ods_SCR[0], raw_ods_TDDR[0], export = export, filename = 'SCR/Subject_1')

plot = False

Subject 2

if plot:
    compare_original_filtered(raw_ods_SCR[1], raw_ods_TDDR[1], export = export, filename = 'SCR/Subject_2')

Subject 3 - Not working

# if plot:
#     compare_original_filtered(raw_ods_SCR[2], raw_ods_TDDR[2], export = export, filename = 'SCR/Subject_3')

Subject 4

if plot:
    compare_original_filtered(raw_ods_SCR[3], raw_ods_TDDR[3], export = export, filename = 'SCR/Subject_4')

Detailed view

%matplotlib inline
sub = 3 # subject number - 1
chan = 20 # channel pick
chan_short = 28 # short channel pick (closest to long channel and same wavelength)
print(raw_ods_TDDR[sub].ch_names[chan])
print(raw_ods_TDDR[sub].ch_names[chan_short])
tmin, tmax = 900, 1000
data_no_SCR = raw_ods_TDDR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SCR = raw_ods_SCR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SC = raw_ods_TDDR[sub].copy().pick(chan_short).crop(tmin,tmax).get_data()[0]
t = np.arange(len(data_no_SCR))/raw_ods[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_SCR, label = 'without SCR', color = 'r')
plt.plot(t, data_SCR, label = 'with SCR', color = 'g')
plt.plot(t, data_SC, label = 'Nearest short channel', color = 'b')
plt.xlabel('time (s)')
plt.ylabel('OD (V)')
plt.title('Detailed view of OD channel S10-D1 760 and short channel S13-D3 760 of subject 4')
plt.legend()
if export:
    fig.savefig('SCR_N_back_subject_4.png')

S10_D1 760
S13_D3 760

%matplotlib inline
sub = 3 # subject number - 1
chan = 8 # channel pick
chan_short = 28 # short channel pick (closest to long channel and same wavelength)
print(raw_ods_TDDR[sub].ch_names[chan])
print(raw_ods_TDDR[sub].ch_names[chan_short])
tmin, tmax = 900, 1000
data_no_SCR = raw_ods_TDDR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SCR = raw_ods_SCR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SC = raw_ods_TDDR[sub].copy().pick(chan_short).crop(tmin,tmax).get_data()[0]
t = np.arange(len(data_no_SCR))/raw_ods[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_SCR, label = 'without SCR', color = 'r')
plt.plot(t, data_SCR, label = 'with SCR', color = 'g')
plt.plot(t, data_SC, label = 'Nearest short channel', color = 'b')
plt.xlabel('time (s)')
plt.ylabel('OD (V)')
plt.title('Detailed view of OD channel S4-D3 760 and short channel S13-D3 760 of subject 4')
plt.legend()
if export:
    fig.savefig('SCR_N_back_subject_4_2.png')

S4_D3 760
S13_D3 760

---

Remark:

• Baseline drifts are still present
• Large oscillation due to physiological noise still present
• Heart beat still present
• SCR did not work as well as for the FTT2

---

Global component derived from the mean or median

It is worth noting that the short channel regression techniques discussed earlier can also be adapted by using the global component derived from the mean or median of all channels. This approach estimates systemic physiological changes common across channels and filters them from each channel, further improving the isolation of neural signals.

Method of Saager and Berger

# Based on code of mne_nirs.signal_enhancement._short_channel_correction
def M_regression_SB(raw, min_dist = 0.01, max_dist=0.045, method = 'mean'):
    raw = raw.copy()

    picks_od = mne.pick_types(raw.info, fnirs='fnirs_od')

    if len(picks_od) == 0:
        raise RuntimeError('Data must be optical density.')
    
    distances = source_detector_distances(raw.info)
    picks_long = picks_od[distances[picks_od] > min_dist]

    Y_long = get_long_channels(raw, min_dist=min_dist, max_dist=max_dist).get_data()
    m = np.zeros((2, Y_long.shape[1]))
    
    if method == 'mean':
        m[0] = np.mean(Y_long[::2], axis = 0) # 760 wavelength
        m[1] = np.mean(Y_long[1::2], axis = 0) # 850 wavelength
    elif method == 'median':
        m[0] = np.median(Y_long[::2], axis = 0) # 760 wavelength
        m[1] = np.median(Y_long[1::2], axis = 0) # 850 wavelength


    for pick in picks_long: # Short channels remain the same!
        i = pick%2
        A_l = raw.get_data(pick).ravel()
        A_s = m[i].ravel()

        # Eqn 27 Scholkmann et al 2014
        alfa = np.dot(A_s, A_l) / np.dot(A_s, A_s)

        # Eqn 26 Scholkmann et al 2014
        raw._data[pick] = A_l - alfa * A_s

    return raw

Estimate global component derived from the mean

# This approach estimates systemic physiological changes common across channels and filters them from each channel
raw_ods_MeanCR = []
for sub in range(4):
    raw_ods_MeanCR.append(M_regression_SB(raw_ods_TDDR[sub].copy(), min_dist=min_dist, method = 'mean'))

Visualize result

Subject 1

if plot:
    compare_original_filtered(raw_ods_MeanCR[0], raw_ods_TDDR[0], export = export, filename = 'MeanCR/Subject_1')

Subject 2

if plot:
    compare_original_filtered(raw_ods_MeanCR[1], raw_ods_TDDR[1], export = export, filename = 'MeanCR/Subject_2')

# Show physiological noise suppresion used for subject 3
plot = True

Subject 3

if plot:
    compare_original_filtered(raw_ods_MeanCR[2], raw_ods_TDDR[2], export = export, filename = 'MeanCR/Subject_3')

plot = False

Subject 4

if plot:
    compare_original_filtered(raw_ods_MeanCR[3], raw_ods_TDDR[3], export = export, filename = 'MeanCR/Subject_4')

Detailed view

%matplotlib inline
sub = 3 # subject number - 1
chan = 20 # channel pick
chan_short = 28 # short channel pick (closest to long channel and same wavelength)
print(raw_ods_TDDR[sub].ch_names[chan])
print(raw_ods_TDDR[sub].ch_names[chan_short])
tmin, tmax = 900, 1000
wavelength = raw_ods_TDDR[sub].ch_names[chan][-3:]
type = 0 if wavelength == 760 else 1
data_no_SCR = raw_ods_TDDR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SCR = raw_ods_MeanCR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_mean = np.mean(raw_ods_TDDR[sub].copy().crop(tmin,tmax).get_data()[type:-4:2], axis = 0)
t = np.arange(len(data_no_SCR))/raw_ods[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_SCR, label = 'without Mean CR', color = 'r')
plt.plot(t, data_SCR, label = 'with Mean CR', color = 'g')
plt.plot(t, data_mean, label = 'mean of all long 760 channels', color = 'b')
plt.xlabel('time (s)')
plt.ylabel('OD (V)')
plt.title('Detailed view of OD channel S10-D1 760 and mean of all long 760 channels of subject 4')
plt.legend()
if export:
    fig.savefig('MeanCR_N_back_subject_4.png')

S10_D1 760
S13_D3 760

%matplotlib inline
sub = 3 # subject number - 1
chan = 8 # channel pick
chan_short = 28 # short channel pick (closest to long channel and same wavelength)
print(raw_ods_TDDR[sub].ch_names[chan])
print(raw_ods_TDDR[sub].ch_names[chan_short])
tmin, tmax = 900, 1000
wavelength = raw_ods_TDDR[sub].ch_names[chan][-3:]
type = 0 if wavelength == 760 else 1
data_no_SCR = raw_ods_TDDR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SCR = raw_ods_MeanCR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_mean = np.mean(raw_ods_TDDR[sub].copy().crop(tmin,tmax).get_data()[type:-4:2], axis = 0)
t = np.arange(len(data_no_SCR))/raw_ods[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_SCR, label = 'without Mean CR', color = 'r')
plt.plot(t, data_SCR, label = 'with Mean CR', color = 'g')
plt.plot(t, data_mean, label = 'mean of all long 760 channels', color = 'b')
plt.xlabel('time (s)')
plt.ylabel('OD (V)')
plt.title('Detailed view of OD channel S4-D3 760 and mean of all long 760 channels of subject 4')
plt.legend()
if export:
    fig.savefig('MeanCR_N_back_subject_4_2.png')

S4_D3 760
S13_D3 760

Estimate global component derived from the median

# This approach estimates systemic physiological changes common across channels and filters them from each channel
raw_ods_MedianCR = []
for sub in range(4):
    raw_ods_MedianCR.append(M_regression_SB(raw_ods_TDDR[sub].copy(), min_dist=min_dist, method = 'median'))

Visualize results

Subject 1

if plot:
    compare_original_filtered(raw_ods_MedianCR[0], raw_ods_TDDR[0], export = export, filename = 'MedianCR/Subject_1')

Subject 2

if plot:
    compare_original_filtered(raw_ods_MedianCR[1], raw_ods_TDDR[1], export = export, filename = 'MedianCR/Subject_2')

Subject 3

if plot:
    compare_original_filtered(raw_ods_MedianCR[2], raw_ods_TDDR[2], export = export, filename = 'MedianCR/Subject_3')

Subject 4

if plot:
    compare_original_filtered(raw_ods_MedianCR[3], raw_ods_TDDR[3], export = export, filename = 'MedianCR/Subject_4')

%matplotlib inline
sub = 3 # subject number - 1
chan = 20 # channel pick
chan_short = 28 # short channel pick (closest to long channel and same wavelength)
print(raw_ods_TDDR[sub].ch_names[chan])
print(raw_ods_TDDR[sub].ch_names[chan_short])
tmin, tmax = 900, 1000
wavelength = raw_ods_TDDR[sub].ch_names[chan][-3:]
type = 0 if wavelength == 760 else 1
data_no_SCR = raw_ods_TDDR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SCR = raw_ods_MeanCR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_median = np.median(raw_ods_TDDR[sub].copy().crop(tmin,tmax).get_data()[type:-4:2], axis = 0)
t = np.arange(len(data_no_SCR))/raw_ods[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_SCR, label = 'without Median CR', color = 'r')
plt.plot(t, data_SCR, label = 'with Median CR', color = 'g')
plt.plot(t, data_median, label = 'median of all long 760 channels', color = 'b')
plt.xlabel('time (s)')
plt.ylabel('OD (V)')
plt.title('Detailed view of OD channel S10-D1 760 and median of all long 760 channels of subject 4')
plt.legend()
if export:
    fig.savefig('MedianCR_N_back_subject_4.png')

S10_D1 760
S13_D3 760

%matplotlib inline
sub = 3 # subject number - 1
chan = 8 # channel pick
chan_short = 28 # short channel pick (closest to long channel and same wavelength)
print(raw_ods_TDDR[sub].ch_names[chan])
print(raw_ods_TDDR[sub].ch_names[chan_short])
tmin, tmax = 900, 1000
wavelength = raw_ods_TDDR[sub].ch_names[chan][-3:]
type = 0 if wavelength == 760 else 1
data_no_SCR = raw_ods_TDDR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_SCR = raw_ods_MeanCR[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_median = np.median(raw_ods_TDDR[sub].copy().crop(tmin,tmax).get_data()[type:-4:2], axis = 0)
t = np.arange(len(data_no_SCR))/raw_ods[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_SCR, label = 'without Median CR', color = 'r')
plt.plot(t, data_SCR, label = 'with Median CR', color = 'g')
plt.plot(t, data_mean, label = 'median of all long 760 channels', color = 'b')
plt.xlabel('time (s)')
plt.ylabel('OD (V)')
plt.title('Detailed view of OD channel S4-D3 760 and median of all long 760 channels of subject 4')
plt.legend()
if export:
    fig.savefig('MedianCR_N_back_subject_4_2.png')

S4_D3 760
S13_D3 760

---

Remark:

• Both techniques yield results similar to those obtained using the SCR.

---

Signal processing methods

Lowpass filtering

In this experiment, mental fatigue is induced through the 2-back task, performed by participants in three consecutive 20-minute blocks (see event ID plot). Over this 1-hour working memory task, we anticipate a significant increase in brain activity in the PFC as participants transition from an alert to a fatigued state (Zhang et al., 2017). The stimulation protocol anticipates sustained brain activity for periods exceeding 100 seconds, suggesting that a bandpass filter with $f_{c, \text{low}} < 0.01$ Hz should be used (Pinti et al., 2019). Consequently, we chose a lowpass filter with $f_{c, \text{high}} = 0.09$ Hz instead of a bandpass filter.

Convert OD data to haemodynamic data

By solving the Beer-Lambert law for two measurement wavelengths on either side of the isosbestic point - the point where the HbR and HbO spectra cross, the unknows can be eliminated from the equation. Small changes of attenuation for both wavelengths, due to changes in HbR and HbO concentrations, can be detected with the following formula: $$ \Delta OD_\lambda = \log \left(\frac{I_{rest}}{I_{test}}\right) \ \approx \epsilon_\lambda^{HbR} \cdot \Delta c^{HbR} \cdot L + \epsilon_\lambda^{HbO} \cdot \Delta c^{HbO} \cdot L $$ From the two resulting equations, one for each wavelength, the chromophore concentrations can be found: $$ [HbO] = \frac{a_{HbO}^{\lambda_2} \cdot \Delta A_{\lambda_1} - a_{HbR}^{\lambda_1} \cdot \Delta A_{\lambda_2}}{L \cdot (a_{HbO}^{\lambda_1} \cdot a_{HbR}^{\lambda_2} - a_{HbO}^{\lambda_2} \cdot a_{HbR}^{\lambda_1})} $$ $$ [HbR] = \frac{a_{HbO}^{\lambda_1} \cdot \Delta A_{\lambda_2} - a_{HbR}^{\lambda_2} \cdot \Delta A_{\lambda_1}}{L \cdot (a_{HbO}^{\lambda_1} \cdot a_{HbR}^{\lambda_2} - a_{HbO}^{\lambda_2} \cdot a_{HbR}^{\lambda_1})}$$

# Convert OD to haemoglobin concentration with build-in function of mne_nirs
# Constant ppf value: PPF = DPF/PVC
raw_haemos_short = []
raw_haemos = []
for sub in range(4):
    if sub != 2:
        raw_haemo = beer_lambert_law(raw_ods_SCR[sub], ppf=6) # ppf=6 is more inline with the community expectations and makes our results easier compared to the majority of the existing literature.
        raw_haemos_short.append(get_short_channels(raw_haemo, max_dist=min_dist)) 
        raw_haemos.append(get_long_channels(raw_haemo, min_dist=min_dist, max_dist=max_dist))
    else:
        raw_haemo = beer_lambert_law(raw_ods_MeanCR[sub], ppf=6) # ppf=6 is more inline with the community expectations and makes our results easier compared to the majority of the existing literature.
        raw_haemos_short.append(get_short_channels(raw_haemo, max_dist=min_dist)) 
        raw_haemos.append(get_long_channels(raw_haemo, min_dist=min_dist, max_dist=max_dist))

Plot the PSD of raw measured data

%matplotlib inline
for sub in range(4):
    fig = raw_haemos[sub].compute_psd().plot(average=True, exclude="bads", xscale='log')
    fig.suptitle('PSD before filtering: subject ' + str(sub+1), weight='bold', size='x-large');
    if export:
        fig.savefig('PSD/PSD_before_filtering_subject_' + str(sub+1) + '.png')

    plt.show()

Effective window size : 204.800 (s)
Plotting power spectral density (dB=True).

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\274252886.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].compute_psd().plot(average=True, exclude="bads", xscale='log')
C:\Users\fabia\AppData\Roaming\Python\Python311\site-packages\mne\viz\utils.py:167: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  (fig or plt).show(**kwargs)

Effective window size : 231.424 (s)
Plotting power spectral density (dB=True).

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\274252886.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].compute_psd().plot(average=True, exclude="bads", xscale='log')
C:\Users\fabia\AppData\Roaming\Python\Python311\site-packages\mne\viz\utils.py:167: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  (fig or plt).show(**kwargs)

Effective window size : 229.376 (s)
Plotting power spectral density (dB=True).

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\274252886.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].compute_psd().plot(average=True, exclude="bads", xscale='log')
C:\Users\fabia\AppData\Roaming\Python\Python311\site-packages\mne\viz\utils.py:167: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  (fig or plt).show(**kwargs)

Effective window size : 251.904 (s)
Plotting power spectral density (dB=True).

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\274252886.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].compute_psd().plot(average=True, exclude="bads", xscale='log')
C:\Users\fabia\AppData\Roaming\Python\Python311\site-packages\mne\viz\utils.py:167: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  (fig or plt).show(**kwargs)

Plot filter response

# Specifications
order = 1000
sample_rates = np.array(sfreqs)  # Example sample rate in Hz
nyquist_rates = sample_rates / 2.0
cutoff_freq = 0.09  # Cut-off frequencies in Hz

# Normalize the cut-off frequencies to Nyquist frequency
normalized_cutoffs = cutoff_freq  / nyquist_rates 

# Design the FIR bandpass filter using firwin
fir_coefficients = []
for sub in range(4):
    fir_coefficients.append(firwin(order + 1, normalized_cutoffs[sub], pass_zero='lowpass'))

# Plot the frequency response
for sub in range(4):
    w, h = freqz(fir_coefficients[sub], worN=8000)
    plt.figure(figsize=(12,5))
    plt.plot(0.5 * sample_rates[sub] * w / np.pi, np.abs(h), 'b')
    plt.title('FIR Lowpass Filter Frequency Response - Subject ' + str(sub+1))
    plt.xlabel('Frequency [Hz]')
    #plt.xscale('log')
    plt.ylabel('Gain')
    plt.xlim((0,0.1))
    plt.vlines(cutoff_freq, -0.1, 1.1, linestyle = '--', color = 'r',  label = 'Cut-off frequencies')
    plt.legend()
    if export:
        plt.savefig('FIR_LPF_response_N_back_subject_' + str(sub+1) + '.png')
    plt.show()

Plot filtered signal

def BPF(raw_input, fir_coefficients, plot = False, export = False, filename = None):

    raw_ = raw_input.copy()
    data = raw_.get_data()
    picks = np.sort(_validate_nirs_info(raw_.info))

    for pick in picks:
        signal = data[pick]
        filtered_signal = filtfilt(fir_coefficients, 1.0, signal) # Zero phase filtering (Can only be done offline)
        
        raw_._data[pick] = filtered_signal
    
    if plot:
        compare_original_filtered(raw_, raw_input, export = export, filename = filename)
    
    return raw_

Subject 1

plot = True

raw_haemos_filtered = []
raw_haemos_filtered.append(BPF(raw_haemos[0], fir_coefficients[0], plot = plot, export = export, filename = 'LPF/N_back_LPF_subject_1'))

plot =False

Subject 2

raw_haemos_filtered.append(BPF(raw_haemos[1], fir_coefficients[1], plot = plot, export = export, filename = 'LPF/N_back_LPF_subject_2'))

Subject 3

raw_haemos_filtered.append(BPF(raw_haemos[2], fir_coefficients[2], plot = plot, export = export, filename = 'LPF/N_back_LPF_subject_3'))

Subject 4

raw_haemos_filtered.append(BPF(raw_haemos[3], fir_coefficients[3], plot = plot, export = export, filename = 'LPF/N_back_LPF_subject_4'))

# Detailed view
%matplotlib inline
sub = 1 # subject number - 1
chan = 4 # channel pick
print(raw_haemos[sub].ch_names[chan])
tmin, tmax = 800, 900
data_no_BPF = 1e6*raw_haemos[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
data_BPF = 1e6*raw_haemos_filtered[sub].copy().pick(chan).crop(tmin,tmax).get_data()[0]
t = np.arange(len(data_BPF))/raw_haemos[sub].info['sfreq'] + tmin
fig = plt.figure(figsize = (12,5))
plt.plot(t, data_no_BPF, label = 'original', color = 'r')
plt.plot(t, data_BPF, label = 'filtered', color = 'g')
plt.xlabel('time (s)')
plt.ylabel('HbO (µM)')
plt.title('Detailed view of HbO channel S3-D1  of subject 2')
plt.legend()
if export:
    fig.savefig('N_back_Detail_BPF_subject_2.png')

S3_D1 hbo

# Possible transient delay if not implementing zero-phase filtering
transient_length = order//2
transient_time = [transient_length/sample_rates[sub] for sub in range(4)]
print('Due to the high order of the filter there is possibly a large transient part of ' + str(transient_time) + ' seconds [' + str(transient_length) + ' samples]')

Due to the high order of the filter there is possibly a large transient part of [50.0, 56.5, 56.0, 61.49999999999999] seconds [500 samples]

for sub in range(4):
    fig = raw_haemos[sub].plot_psd(average=True, fmax=2, xscale='log', color='k', show=False);
    raw_haemos_filtered[sub].plot_psd(average=True, fmax=2, xscale='log', ax=fig.axes, color='g', show=False);
    leg_lines = [line for line in fig.axes[0].lines if line.get_linestyle() == '-']
    fig.legend(leg_lines, ['Filtered Data', 'Measured Data'], loc="lower left", bbox_to_anchor=(0.15, 0.2))
    fig.suptitle('Comparison PSD of filtered and unfiltered haemodynamic data of subject ' + str(sub + 1), fontsize=16, fontweight='bold')
    if export:
        fig.savefig('PSD/N_back_PSD_comparison_after_BPF_subject_' + str(sub+1) + '.png')

NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 204.800 (s)
Plotting power spectral density (dB=True).
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 204.800 (s)
Plotting power spectral density (dB=True).
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 231.424 (s)
Plotting power spectral density (dB=True).

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:2: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].plot_psd(average=True, fmax=2, xscale='log', color='k', show=False);
C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  raw_haemos_filtered[sub].plot_psd(average=True, fmax=2, xscale='log', ax=fig.axes, color='g', show=False);
C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:2: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].plot_psd(average=True, fmax=2, xscale='log', color='k', show=False);

NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 231.424 (s)
Plotting power spectral density (dB=True).
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 229.376 (s)
Plotting power spectral density (dB=True).
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 229.376 (s)
Plotting power spectral density (dB=True).
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().
Effective window size : 251.904 (s)
Plotting power spectral density (dB=True).
NOTE: plot_psd() is a legacy function. New code should use .compute_psd().plot().

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  raw_haemos_filtered[sub].plot_psd(average=True, fmax=2, xscale='log', ax=fig.axes, color='g', show=False);
C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:2: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].plot_psd(average=True, fmax=2, xscale='log', color='k', show=False);
C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  raw_haemos_filtered[sub].plot_psd(average=True, fmax=2, xscale='log', ax=fig.axes, color='g', show=False);
C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:2: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  fig = raw_haemos[sub].plot_psd(average=True, fmax=2, xscale='log', color='k', show=False);

Effective window size : 251.904 (s)
Plotting power spectral density (dB=True).

C:\Users\fabia\AppData\Local\Temp\ipykernel_15180\3276914840.py:3: FutureWarning: The value of `amplitude='auto'` will be removed in MNE 1.8.0, and the new default will be `amplitude=False`.
  raw_haemos_filtered[sub].plot_psd(average=True, fmax=2, xscale='log', ax=fig.axes, color='g', show=False);

---

Remark:

• High frequency content of the signals removed
• Use zero-phase filtering as we are working with offline signals

---